Wave packet dynamics and valley filter in strained graphene 
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The time evolution of a wavepacket in strained graphene is studied within the tight-binding model 
and continuum model. The effect of an external magnetic field, as well as a strain-induced pseudo- 
magnetic field, on the wave packet trajectories and zitterbewegung are analyzed. Combining the 
effects of strain with those of an external magnetic field produces an effective magnetic field which 
is large in one of the Dirac cones, but can be practically zero in the other. We construct an efficient 
valley filter, where for a propagating incoming wave packet consisting of momenta around the K 
and K' Dirac points, the outgoing wave packet exhibits momenta in only one of these Dirac points, 
while the components of the packet that belong to the other Dirac point are refiected due to the 
Lorentz force. We also found that the zitterbewegung is permanent in time in the presence of either 
external or strain-induced magnetic fields, but when both the external and strain-induced magnetic 
fields are present, the zitterbewegung is transient in one of the Dirac cones, whereas in the other 
cone the wave packet exhibits permanent spatial oscillations. 

PACS numbers: 72.80.Vp,73.23.-b, 85.75.Hh 



I. INTRODUCTION 

Since its first synthesis in 2004fi graphene has been at- 
tracting much interest due to its unique electronic prop- 
erties arising from its singular energy spectrum, where in 
the vicinity of the points labelled as K and K' in recip- 
rocal space, the charge carriers behave as massless quasi- 
particles and exhibit an almost linear dispersion.- These 
quasi-particles obey the Dirac- Weyl equations and there- 
fore should be subject to quasi-relativistic effects, such as 
zitterbewegung, i.e., a trembling motion caused by inter- 
ference between positive and negative energy states<^"— 
The phenomenon of zitterbewegung was predicted in 
1930 by Schrodinger— and has been subject of renewed 
interest over the past years. Previous theoretical works 
have suggested few ways of experimentally observing 
zitterbewegung, e. g. in narrow-gap semiconductors^, 
in III-V zinc-blende semiconductor quantum wells with 
spin-orbit coupling^ and, more recently, in monolayer^ 
and bilayer^'^ graphene. An experimental simulation of 
the zitterbewegung of free relativistic electrons in vacuum 
was performed by Gerritsma et a\i^ by using trapped 
ions. 

Strain engineering in graphene has recently become 
a widely studied topicJ^"— The elastic properties of 
graphene nanoribbons were theoretically investigated by 
Cadelano et alM-, which studied the in-plane stretch- 
ing and out of plane bending deformations by combining 
continuum elasticity theory and tight-binding atomistic 
simulations. Later, Cocco et alJ^ and Lu and Guc^" 
showed that a combination of shear and uniaxial strain 
at moderate absolute deformations can be used to open 
a gap in the graphene energy spectrum. It has been 
shown recently that specific forms of strain produce a 
pseudo-magnetic field in graphene, which does not break 



the time reversal symmetry and which points in opposite 
directions for electrons moving around the K and K' 
pointsr ^ The strain-induced magnetic field is expected 
to produce Landau levels and, consequently, the quan- 
tum Hall effect, even in the absence of an external mag- 
netic field/S^'SS, Guinea et al.^ showed theoretically that 
an in-plane bending of the graphene sheet leads to an 
almost uniform field. Landau levels as a consequence of 
strain-induced pseudo-magnetic fields greater than 300 
Tesla were recently observed with scanning tunnelling 
microscopy (STM) in nanometer size nanobubblesj^ 

Although previous works have studied wave packet 
propagation for the Dirac- Weyl Hamiltonian of graphene 
in the absence of external fields and potentials)^ or 
in the presence of uniformly applied external magnetic 
fields^i^i"— , there is still a lack of theoretical works on the 
wave packet propagation through potential and (pseudo) 
magnetic field step barriers. Moreover, the time evolu- 
tion of a wave packet in graphene within the tight- binding 
(TB) model, where the intra-valley scattering to higher 
energy states and inter-valley scattering due to defects 
appear naturally, is still hardly found in the literature. 
It is also interesting to see whether the results from Dirac 
and TB approaches for graphene differ or are similar. 

In this paper, we investigate the time-evolution of 
wavepackets in graphene within the TB model, based 
on the split-operator technique for the expansion of the 
time-evolution operator. We trace a parallel between 
the results from the TB model and those obtained from 
the Dirac approximation for particles with momentum 
close to one of the Dirac cones of the Brillouin zone 
of graphene. The proposed method is then applied to 
the study of the dynamics of Gaussian wave packets in 
graphene under external magnetic fields. The effects 
of the pseudo-magnetic field induced by bending the 
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graphene sheet into an arc of circle are analyzed as well. 
Our results show that for an appropriate choice of strain 
and external magnetic field strength, the system exhibits 
a strong effective magnetic field for particles in one of the 
Dirac cones, whereas in the other cone the external and 
pseudo-magnetic fields cancel each other and the effective 
magnetic field is practically zero. We show that this effect 
is manifested as a transient (permanent) zitterbewegung 
for electrons in the cone where the effective magnetic field 
is zero (non-zero), which can be verified experimentally 
by detecting the electric field radiation emitted by the 
trembling wave packet^l Moreover, our results show that 
with such a choice of external and strain-induced mag- 
netic fields, one can construct an efficient valley-filter, 
which can be useful for future valley-tronic devices.— 



II. TIME EVOLUTION OPERATOR 

By solving the time-dependent Schrodinger equation, 
one obtains that the propagated wavefunction after a 
time step At can be calculated by applying the time- 
evolution operator on the wave packet at any instant t 



*(f,t + At) = e" 



rHAt 



^'(f,i). 



(1) 



where we assumed that the Hamiltonian H is time- 
independent. Different techniques have been suggested 
for the expansion of the exponential in Eq. ([T]), e.g. the 
Chebyschev polynomials method^ and the second order 
differencing scheme^^. The numerical method that we 
use for the application of the time evolution operator in 
this work, namely, the split-operator method,— is the 
subject of this section, where we will show how this tech- 
nique can be adapted for the study of the wave packet 
dynamics in graphene within the tight-binding and Dirac 
approximations. 



A. Tight-Binding model 

Graphene consists of a layer of carbon atoms form- 
ing a honeycomb lattice, which can be described by the 
Hamiltonian 



Htb — 



E 

<i,j> 



TijCiCj 



(2) 



where Ci(cj) annihilates (creates) an electron at the site i, 
with on-site energy e^, and the sum is taken only between 
nearest neighbor sites i and j, with hopping energy Ty. 
The effect of an external magnetic field can be calculated 
by including a phase in the hopping parameters accord- 
ing to the Peierls substitution Tij — >■ Tij exp 



h J J 



dl 



where A is the vector potential describing the mag- 
netic fieldj^"^ In a strained graphene sheet, the dis- 
tance between two adjacent sites i and j is changed by 



unstrained graphene and aij is the distance between the 
sites after the stress. The change in the inter-sites dis- 
tance affects the hopping energy between the sites, which 
becomes^! Tij — >■ Ty (1 -I- 2Aajj/ao). A similar expres- 
sion can be obtained by expanding Eq. (13) of Ref. 
[37I in Taylor series and neglecting higher order terms in 
Aoij, i.e., considering small lattice deformations. The 
strain-induced change in the hopping energies leads to 
an effective pseudo-magnetic field, which points to op- 
posite directions in the valley K and K' ?^ Notice that 
the pseudo-magnetic field in our model is not introduced 
artificially by considering an additional vector potential 
in the Peierls phase, but it rather appears naturally after 
the changes in the inter-site distances due to the strain. 



a) {l,l}^{i,2} 




FIG. 1: (Color online) (a) Sketch of the honeycomb lattice of 
graphene, made out of two superimposed triangular lattices 
A and B. The atoms are labelled as {n, m} according to 
their line and column numbers n and m, respectively, (b) 
Reciprocal lattice of graphene, with K (black) and K' (gray) 
Dirac points, where the area defined by the reciprocal vectors 
61 — (— 27r/3ao, \/3) and &2 ~ (47r/3ao,0) represents the first 
Brillouin zone. The numbers close to each Dirac point are 
explained in the text. 

Let us label the sites i of the graphene lattice according 
to their line and column numbers n and to, respectively, 
as shown in Fig. [IJa). The basis vector state |n, to) 
represents an electron confined on the site of line n and 
column TO. In a lattice consisting only of non-interacting 
sites, each \n^m) is an eigenstate of Htb with energy 
En.m, i-e. HTB\n,m) — en,m\n,m) . Limiting ourselves to 
nearest neighbors interactions, we find 

HTB\n,m) = enm\n,m) +T,n+i\n,m + 1) 
+r„i_i|n,TO - 1) + Tn+i\n + 1,771) + T„_i|7i - 1,to),(3) 

where T„±i and Tm±i are equivalent to the hopping en- 
ergies Tij between the site i = {n, to} and the adjacent 
sites j = {77± 1, to} and j — {n, to± 1}, respectively. Eq. 
^ can be rewritten as 



HTB\n, ni) = Hn\n, m) + i?m|r7, m), 
where the operators Hn and _ff,„ are defined as 



(4) 



Aoi 



ao, where ag is the lattice parameter of 



Hn\n,m) =Tm+i\n,m+l)+Tm-i\n,m-l) + ^Y- (5a) 
and 

Hjn\n,m) = Tn+i\n+l,m)+Tn-i\n-l,m) + ^^. (5b) 
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The wavefunction at any instant t is then written as a 
hnear combination of the basis vector states ^'^ ^ = 
X)n m ^n.ml"' "^)- ^hc advantage of following the pro- 
cedure described by Eqs. (|3][S|) lies in the fact that the 
operators Hn and Hm in Eq. ([5]) can be represented by 
tri-diagonal matrices, which are easier to handle than the 
matrix representing the full Hamiltonian, Eq. (|3]). 

The split-operator technique can now be applied to the 
Hamiltonian Eq. (jlj , so that the time evolution operator 
is approximated by 



(6) 

where the error comes from the non-commutativity be- 
tween the operators 77„ and Hm- We drop the 0{At^) 
terms by considering a small time step At = 0.1 fs. The 
propagated wavefunction is then obtained from Eq. ([T}, 
which in this case reads 

^t+At ^ ^^H„,At -^j-H„At -^H„,At^t (j-. 
^ n.rn c- o ^ n,m- \'J 

This equation is solved in three steps: 

Un.m. — e -"^ 



'rH„At^ 



t+At _ -T^H^At 



(8a) 



(8b) 



(8c) 



Using the Cayley form for the exponentials,— we can 
rewrite Eq. (|S^) as 



-H^At.j.t 



1 iAt TT 

^ 4R ,T,t 



1 + 



vf* +0(At^),(9) 



which leads to 



1 TT 1 

1 + ] Tin 



(10) 



As the wavefunction Vt^ ^ is already known, the matrix 
equation in Eq. (ITUl) can be straightforwardly solved to 
obtain r]n,m- We repeat this procedure for the other two 
exponentials in Eqs. and (Ht), and eventually obtain 

xUt+At 

In fact, the form in Eq. ([TU]) can also be applied to the 
full Hamiltonian Htb, i-e. without splitting the Hm and 
Hn terms. However, this would lead to matrix equations 
involving penta-diagonal matrices, which are harder to 
handle than the tri-diagonal matrices in Eq. PH)) . As the 
error produced by the splitting in Eq. ([7]) is smaller than 
the error produced by the (necessary) expansion of the 
exponential given by Eq. (jlOp , it is worthy to split these 
terms in order to simplify the numerical calculations. 



B. Dirac-Weyl equation 

From the TB Hamiltonian Eq. ([2]), considering an in- 
finite graphene sheet in the absence of external potential 
and magnetic fields, one obtains the energy bandstruc- 
ture of graphene as 



E{k] =±TJ3 + f{k 



0{At^), f (k^ = 2 cos (V3fcyao)-l-4cos 



— kyao I cos I -k^ao 



(11) 

which is gapless in six points of the reciprocal space where 
E = 0, from which only two are inequivalent, labelled as 
K and K', as shown in Fig. [IJb)j^i2S In the vicinity 
of each of these points, the dependence of the energy 
spectrum on the wave vector k is almost linear and the 
electron can be described as a massless fermion by the 
Dirac Hamiltonian 



Hd = 



vpS ■ (p + eA) + y(x,y)I e 



-20 



(12) 



where vp = 3Ta()/2h is the Fermi velocity, A is the elec- 
tromagnetic vector potential, V{x,y) is an external po- 
tential, I is the identity matrix, a is the Pauli vector and 
the wavefunctions are pseudo-spinors ^E* = [^'aj^b]"^, 
with 'i>A(B) the probability of finding the electron in 
the sub-lattice A{B) that composes the honeycomb lat- 
tice of graphene*^ The angle </> is different for electrons 
around the K and K' Dirac cones. In the vicinity of the 
k-th Dirac point (see labels for each Dirac point in Fig. 
[IJb)), one obtains (j) = — 7r/6 + kw/S, with k = 1 — 6. 
From here onwards, we will refer to the coordinates in 
the Dirac (TB) model as x (x) and y (y). 

The exponential term in Eq. is usually dropped, 
because it can be considered as a phase in the state vec- 
tors in the Dirac model. However, this term has an im- 
portant meaning when comparing with the TB model: 
this exponential can be identified as a rotation operator 
with angle (jj- Notice that an infinite graphene hexago- 
nal lattice has C^y symmetry, i.e. it is symmetric only 
for rotation angles fc7r/3 (fc - integer). As a consequence, 
the Dirac Hamiltonian Eq. ([T^ without the exponen- 
tial term is not symmetric in the momenta px and py, as 
already pointed out previouslyjS, So, what would be the 
meaning of the direction-dependent observables in the 
Dirac description of graphene, when they are not sym- 
metric under rotation, exhibiting a privileged direction? 
Defining y (x) as the zigzag (armchair) direction, as in 
Fig. [ija), the results obtained by the Dirac approxima- 
tion for the X and y components of any observable are 
compared to the armchair and zigzag directions, respec- 
tively, after performing a rotation (p in the coordinates 
of the Dirac model. From the possible values of cj), one 
deduces straightforwardly that at any Dirac cone, the co- 
ordinate X (y) of the Dirac model is always related to one 
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of the zigzag (armchair) directions of the real graphene 
lattice. On the other hand, for finite rectangular sam- 
ples the different angles 4> represent two distinguishable 
situations, since the rectangle does not share the C^v 
symmetry of the infinite graphene lattice: the x direction 
in the Dirac model for the odd (even) cones in Fig. [ijb) 
represents a diagonal (vertical) direction in the rectangle. 
The comparison between the TB model for a rectangu- 
lar graphene flake and the Dirac approximation will be 
discussed in details further, in the following section. 

In a recent work, Maksimova et al^ presented an an- 
alytical study of the time evolution of a Gaussian wave 
packet in graphene in the absence of external potentials 
and/or magnetic fields within the continuum model, i. 
e. using the Dirac- Weyl Hamiltonian for electrons in the 
vicinity of the Dirac point K. In this paper, we will 
use an alternative and more general way of calculating 
the dynamics of a wavepacket in such a system,^ based 
on the split-operator technique, which can be applied for 
systems under arbitrary external potentials and magnetic 
fields. 

The Dirac- Weyl Hamiltonian Hu in Eq. (fT2|) can be 
separated as -ffu Hk + Hr, where Hk — hvpa ■ k 
keeps only the terms which depend on the wave vector 
k, whereas Hr = vpea ■ A + VI depends only on the real 
space coordinates x and y. Following the split-operator 
method, the time evolution operator for the Hamiltonian 
Hu can be approximated as 



exp 



'-^{Hk+Hr) 













iAt 


exp 




exp 


^fc 

h 


exp 


Hr 

2h 



(13) 



with an error of the order of 0{At^), due to the non- 
commutativity of the operators involved. We invoke a 
well-known property of the Pauli vectors 



exp 



=cos(5)I-z5^(5.<7 



S 



(14) 



for any vector S, where S — \S\, and rewrite the ex- 
ponentials in real and reciprocal space, respectively, in 
matrix form 



M, 



cos (A) 1 — i 



. sin (A) 



Ax 





+ iAv 



Mfc — cos(/t)I — i 



sin(fi;) 



Ax ^ Aw 



e 2h 
(15a) 

(15b) 



where A = |A| = AtvFe\A\/2h, k, = Ato^k and k = 



|k| = Atvpy^y^ + ky, so that the time evolution of a wave 

packet '^i:i{x,y) — [0a, 0s]^5'(a;, y) can be calculated as 
a series of matrix multiplications 



The matrix multiplication by Mfc is made in reciprocal 
space by taking the Fourier transform of the functions. 
In the absence of magnetic fields and external potentials, 
one has = I and 



«'(r,t + At) =Mfe*(r,t), 



(17) 



where the matrix multiplication in reciprocal space gives 
the exact result for the time evolution of the wave packet, 
since there is no error induced by non-commutativity of 
operators or matrices in this case. This shows that the 
split-operator method provides a way to study the dy- 
namics of wavepackets in graphene within the continuum 
model where, in the presence of magnetic fields and/or 
external potentials, one can control the accuracy of the 
results by making At smaller, while in their absence, the 
problem is solved exactly by a simple matrix multiplica- 
tion for any value of At. 



III. RESULTS AND DISCUSSION 



We shall now discuss the results obtained for a 
graphene lattice with 2000x3601 atoms, with armchair 
(zigzag) edges in the a;(2/)-direction. The nearest neigh- 
bors hopping parameter and the lattice constant of 
graphene are taken as r = —2.7 eV and oo — 1.42 A, 
respectively. 

As initial wave packet, we consider a Gaussian centered 
at ro = (a;o, yo) in real space and q = {q^, q^) in reciprocal 
space: 



exp 



(x - xpf + [y - yoY 
2d? 



+ iq ■ f 



4'(r, t + At) ^M.r-M.k- M^*(r, t) + 0{At^ 



(16) 



(18) 

where N is the normalization factor. Notice that we 
have included a pseudo-spinor [ci, in the initial wave 
packet, where Ci(2) is the probability of finding the elec- 
tron in the triangular sub-lattice A{B) that composes 
the graphene hexagonal lattice. One can also rewrite the 
pseudo-spinor as [l,e'^]"^, where the pseudo-spin polar- 
ization angle 6 is shown explicitly. The pseudo-spin is a 
concept normally attributed to the Dirac description of 
graphene. Indeed, the pseudo-spin of a wavefunction in 
the Dirac model is related to the expectation values of 
the Pauli matrices (ci), which can involve integrals of the 
product between wavefunctions for sub- lattices A and B. 
Such a definition fails for the TB wavefunctions, since in 
this case they are defined in different points of the lat- 
tice, so that any integral that mixes functions of both 
sub-lattices gives zero. Even so, the study of the pseudo- 
spin related to the initial discrete wave packet helps to 
understand the wave packet dynamics obtained from the 
TB model, as we will see further in this section. 



5 



A. Initial pseudo-spin polarization and 
zitterbewegung revisited 

In this subsection, we will use the TB model to review 
some of the main properties of the wave packet dynam- 
ics in graphene. Within the TB model, we consider the 
initial wave packet as a discrete form of the Gaussian dis- 
tribution in Eq. ([T8|) for the graphene hexagonal lattice, 
where we multiply the Gaussian function by Ci(2) in the 
sites belonging to the triangular sub-lattice A{B). From 
Eq. dill) , it is clear that in momentum space, the region 
of interest is the vicinity of the Dirac points K and K' , 
since the energy corresponding to wave vectors out of this 
region is very high. 



, 400 
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x(A) 



20 40 60 80 100 
t(fs) 



FIG. 2: (Color online) (a) Contour plots of the squared mod- 
ulus of the wavefunction after a time evolution of t = 40 fs, 
for three different values of the dimensionless parameter kyd. 
(b) Expectation value of the position and (c) velocity in y- 
direction as a function of time. The results obtained from 
the TB (Dirac) model are presented as curves (symbols), for 
kyd = 1 (solid, circles), 2 (dashed, triangles) and 4 (dotted, 
squares) 



In the TB model for two-dimensional crystals, one 
usually considers the same Gaussian distribution for all 
the sites of the lattice.'"^ This is equivalent as choosing 
Ci = C2 = 1 in Eq. (ITSt . Figure [2][a) shows the contour 
plots of the squared modulus of the propagated wave- 
functions at t = 40 fs for these values of Ci, considering 
an initial wave vector q = (0, ky) + K, i.e. in the vicinity 
of the K point labelled as 2 in Fig. [Ijb). As shown in 
Ref. [2^ the wave packet dynamics near the Dirac cones 
in graphene does not depend separately on the momen- 
tum ky or on the width d, but rather on the dimensionless 



quantity kyd. This result was obtained from the Dirac 
model for graphene, i.e. considering that even high en- 
ergies states exhibit linear dispersion. Within the TB 
model we expect that wave packets with the same kyd 
behave alike only if ky is not too far from the Dirac cone 
and if d is not too small, so that the packet is well lo- 
calized in energy space. Within these conditions. Fig. [5] 
shows the time evolution for different values of this di- 
mensionless quantity: kyd — 1 and 2, with d = 100 A , 
and k^d = 4, with d = 200 A . We observe that the dis- 
persion of the wave packet is stronger for smaller values 
of kyd, where it becomes distorted into an arc-like shape. 
For larger kyd, on the other hand, the wave packet keeps 
its circularly symmetric shape for longer times. 

As explained in the previous subsection, in order to 
obtain the Dirac Hamiltonian Eq. one must shift 

the origin of the wave vectors k to one of the six Dirac 
points shown in Fig. [IJb). Besides, one must also rotate 
the X and y axis by an angle (j) which depends on the K 
or K' point that is considered as the origin in momentum 
space. For the K — {0, 47r/3-\/3ao) point, labelled as 2 in 
Fig. [IJb), the Dirac Hamiltonian Eq. (IT^ is obtained by 
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FIG. 3: (Color online) Final velocities for the Gaussian wave 
packet in Eq. (|18[) . with pseudo-spin ci = C2 = 1 and mo- 
mentum q = {0, ky) + K as a function of (a) the momentum 
A:^, for widths d = 50 A , 100 A and 200 A , and (b) the width 
d, for momenta k° = 0.01 A'S 0.02 A and 0.04 A The 
symbols (curves) are obtained from the TB (Dirac) model. 
In (b), the results from the Dirac model for = 0.02 A"^ 
(dashed) and 0.04 A~^ (solid) are multiplied by 0.985 and 
0.97, respectively. 
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rotating the axis by 90°, with other words, by a trans- 
formation of the coordinates as x — >■ — y and ?/ — >■ x. The 
pseudo-spinor ci — C2 — I represents a wave packet po- 
larized in x-direction, i.e. (ax) > and (a^) = (ay) = 0. 
From the Heisenberg picture, we obtain the velocity in 
the x-direction for the proposed wave packet as 



just for completeness: 



dx 
di 



D\ 



VpCTx- 



(19) 



Performing the appropriate coordinate transformations, 
the velocity obtained from the TB model for the y- 
direction must be consistent with the prediction from the 
Dirac approximation, namely, Vy — dx/dt = vpUx. This 
suggests that such a wave packet propagates towards the 
positive y-direction, but with non-constant velocity, since 
ax does not commute with Hjj. The expectation value 
of the y position of the packet (y) is shown as a function 
of time by the curves in Fig. [2][b), for fc°d = 1 (sohd), 2 
(dashed) and 4 (dotted), where the results obtained by 
the Dirac equation are shown as symbols for comparison. 
A different linear behavior is already observed for each 
wave packet at large time, implying that they have dif- 
ferent velocities, which is kind of counter-intuitive, since 
low-energy electrons in graphene are expected to propa- 
gate always with the same Fermi velocity vp- Figurell^c) 
shows the velocity Vy , calculated by taking the derivative 
of the TB results for (y) with respect to time, which ex- 
hibits clear oscillations that are damped as time evolves, 
converging to a final value v^^ < vp that depends on 
the initial wave packet's width d and wave vector ky. 
The velocities obtained by the Dirac model are shown by 
symbols, where the same qualitative behavior is observed 
as obtained from the TB model, though for higher wave 
packet momentum and width, a small quantitative differ- 
ence is observed, which is a consequence of the different 
energy-momentum dispersion. 

The oscillatory behavior of the velocity is a manifesta- 
tion of the zitterbewegung, i.e. a trembling motion of the 
wave packet due to the interference between positive and 
negative energy states that makes up the initial Gaussian 
wave packet This effect is well-known for relativistic 
particles, which are described by the Dirac equation, and 
is also observed for electrons in graphene in the vicinity 
of the K and K' points, since they can be described as 
massless quasi-particles by the Dirac equation as well. 
The velocity wiggles with shorter period and smaller am- 
plitudes for larger values of kyd. The convergence of the 
velocities demonstrates that the zitterbewegung is not a 
permanent, but a transient effect 

Figure [3] shows the converged velocity as a func- 
tion of (a) the momentum ky and (b) the width d of 
the Gaussian wave packet. The TB results (symbols) are 
compared to those calculated from Eq. (31) in Ref. [2^ 
(curves) , which was obtained analytically from the Dirac 
approximation in the i — > cx) limit and is repeated here 



^ = 1 - 
VF 



2{k^Jf 



(20) 



Within the Dirac model, one can observe that increasing 
d or in Eq. (|20p . the final velocity increases mono- 
tonically and approaches vp, which is reasonable, since a 
wider packet in real space leads to a narrower distribu- 
tion in fc-space, whereas a higher value of the wave vector 
makes the center of the packet lay far from E = Q. In 
both cases, the interference with negative energy states 
is reduced and, consequently, zitterbewegung effects are 
less significant. However, since the analytical formula 
Eq. (I^Hl) does not take into account any effect such 
as the curvature of the energy bands for higher energy 
states or trigonal warping effects, this formula is not ex- 
pected to give accurate results for larger ky. Indeed, Fig. 
|31[a) shows that a very good agreement between TB and 
Dirac results can be observed only for small values of 
whereas for large fcJJ, the final velocities obtained from 
the TB model are lower than those obtained from the 
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FIG. 4: (Color online) (a) Contour plots of the squared modu- 
lus of the wavefunction after a time evolution of f = 40 fs, for 
three initial configurations of pseudo-spin [ci,C2]"^ and mo- 



mentum go: 1) [1,0]"^, k'i 



and kyd - 



4; 2) fc!^ 







and kid = 4; and 3) k°d = 4, = 0. (b) Expec- 

tation value of X obtained by the TB model for the initial 
wavepackets 1 (solid), 2 (dashed) and 3 (dotted) as a func- 
tion of time. The results obtained by the Dirac model, after 
the appropriate coordinates rotation (see text), are shown as 
circles, triangles and squares, respectively. The inset shows 
the trajectory of the wavepacket obtained from the TB model 
for the initial wavepacket 3. 
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FIG. 5: (Color online) (a) Expectation value of the velocity 
as a function of time, for wavepackets with ky = — 
and pseudo-spinor [1,1]"^ (solid) and [1,*]"^ (circles) at the 
Dirac point K = (0, 47r/3%/3ao) (point 2 in Fig. [Hb)), and 
[1,1]"^ (triangles) at K' = (27r/3ao, 27r/3\/3ao) (point 1 in 
Fig. [Hb)). The x and y components of the velocity in the 
latter case are shown as dashed and dotted curves, respec- 
tively. 



Dirac model and do not increase monotonically, but de- 
creases slowly for very large fc^ , as a consequence of the 
curvature of the energy bands. On the other hand, in 
Fig. Efb) we observe that varying the wave packet width 
for a fixed momentum, good qualitative agreement with 
the Dirac model is obtained for almost any value of d. 
The curves for larger values of ky (solid and dashed) are 
just quantitatively different from those obtained by the 
TB model, and they are comparable to the TB results af- 
ter multiphcation by a factor 0.985 (0.970) for k° = 0.02 
(0.04 A~^). Worse qualitative agreement between 
TB and Dirac results in this case is observed only for 
very small d, where the Gaussian width in energy space, 
given by AE = vpfidr^, incorporates higher energy val- 
ues, leading to deviations in obtained from the TB 
model as compared to those from the Dirac model. 

In Fig. mja) we show contour plots of the squared 
modulus of the wave function at t = 40 fs for three dif- 
ferent choices of wave vectors q= {k'^,ky)+K and initial 
pseudo-spinors: 1) [1,0]-^, with = and kyd = 4, 2) 
with kl = and k^d = 4, and 3) with 
k^d = 4 and fc" = 0. The curves (symbols) in Fig. 
mjb) show the expectation value {x) for each case ob- 
tained by the TB (Dirac) model. In case 1 (solid, cir- 
cles), the pseudo-spinor points in the z-direction, so that 
(cs) = (cy) and, consequently, the velocity for both in- 
plane directions must be zero. Indeed, the wave packet 
splits into two equal parts that propagate in opposite y 
directions, leading to Vy = 0. In the x-direction, although 



there is still a small zitterbewegung, {x) rapidly converges 
to a constant, leading to — 0. In case 2 (dashed, tri- 
angles), the pseudo-spinor points in the y-direction, but 
the momentum of the wave packet in this direction is 
zero, so that the packet splits in the y-direction, since 
Vy — vpf^x = 0, but drifts slowly in the —x direction (or, 
equivalently, y). In case 3 (dotted, squares), both the 
pseudo-spin and the momentum are in the y-direction, so 
that the wave packet propagates in this direction without 
any splitting. This situation is comparable to the one in 
Fig. I^a), since in both cases the pseudo-spin and mo- 
mentum are in the same direction and, as a consequence, 
the wave packet propagates in this direction practically 
preserving its circular symmetry. However, in the case 
3, the packet still presents a very weak oscillation in the 
y-direction, and also drifts very slowly in this direction, 
as one can see from the trajectory of the packet in the 
X — y plane for this case, shown in the inset of Fig. Hl^b). 
This oscillation and drift are related to the contributions 
of higher energy states in the wave packet: a wave packet 
centered around fci"' = and k'^^ ^ 0, as in Fig. HJa), 
has a symmetric distribution of momenta in x-direction 
even for higher energies and, consequently, there is no ad- 
ditional oscillation in this direction. On the other hand, 
a packet centered around k'^x^ ^ and = 0, as in Fig. 
Sfb) , does not have a symmetric distribution of momenta 
in the y-direction due to the trigonal warping for higher 
energies and, consequently, some oscillations are observed 
in this direction. As the standard Dirac Hamiltonian 
for graphene does not take trigonal warping into account, 
this effect is not observed in the Dirac model. 

In the numerical work of Thaller,^ as well as in the an- 
alytical work of Maksimova,^^ it is demonstrated within 
the Dirac model that even when ky — k^ — 0, wave 
packet motion is still observed due to zitterbewegung ef- 
fects. The velocities of the wave packet obtained from 
our TB model of graphene for wave packet momenta ex- 
actly at the K' and ivT, i.e. points 1 and 2 in Fig. [TJa), 
respectively, are shown in Fig. [5l The velocities exhibit a 
damped oscillation with the same time-dependent mod- 
ulus for any pseudo-spin and Dirac point, though they 
point in different directions: for [1,1]"^ (solid) and [1, i]"^ 
(circles) in K, the velocity is in the y and —x direc- 
tions of the lattice, respectively, which are exactly the 
directions of polarization of these pseudo-spinors after 
the = 7r/2 rotation required by the K cone 2. In the 
K' cone 1, the rotation angle is (f) = 7r/6 and, accord- 
ingly, the velocity points in this direction, as one can see 
by the decomposition of the velocity in the components 
{vx) (dashed) and {vy) (dotted), which obey exactly the 
relations {vx) = (\/3/2)(w) and {vy) = (l/2)(w). Notice 
that the velocities converge exactly to u_f/2, a value that 
can also be obtained analytically by making kyd — > in 

Eq. dani). 

Henceforth, we will consider initial wave vectors (fo 
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FIG. 6: (Color online) (a) Sketch of the strained graphene 
sheet: we consider a rectangular sample of width W and 
height L, bent into an arc of circle with inner radius R. The 
unstrained graphene sheet is shown as open circles, for com- 
parison, (b) Strain-induced magnetic field barrier step, ob- 
tained by bending the graphene lattice only in the y > 
region. The number of atoms was reduced in both figures, in 
comparison to the lattices studied in this work, in order to 
improve the visualization. 

around the Dirac points 2 and 5 of Fig. [Ijb), namely 

K=(o,^^] and K'=(o, (21) 

V 3v3ao/ \ 3v3ao/ 

respectively. This choice is very convenient, since the ro- 
tation angles for these points are (p — 7r/2 and 37r/2, re- 
spectively, so that the pseudo-spinor [1,1]"^ points to the 
y {-y) direction in the former (latter) case. Hence, with 
this pseudo-spinor, wave packets in K (K') will propa- 
gate with positive (negative) velocity in the vertical zig- 
zag direction. 

B. External magnetic fields and strain 

Recently}^ it was shown theoretically that bending a 
graphene sheet into an arc of a circle produces a strong 
and almost uniform pseudo-magnetic field profile. Fig- 
ure [Hlja) illustrates such a strained system, where the 
rectangular graphene sample of width W and height L is 
bent into an arc of a circle with inner radius R. As the 
(pseudo) magnetic field points in the same direction (op- 
posite directions) at each K and K' points, the combi- 
nation of both external and strain-induced magnetic field 
effects provides a valley-dependent magnetic field. If one 
applies the appropriate external magnetic field for some 
configuration of the strained graphene, one can obtain 
an almost perfect suppression of the effective magnetic 
field at one of the Dirac cones, while the effective field 
in the other cone is enhanced. This leads to a compli- 
cated system to be studied within the Dirac approxima- 
tion, since one has two completely different systems for 



the K and K' valleys. Namely, Landau levels would be 
present only around one of the cones (though one can- 
not expect a perfect Landau level spectrum, since the 
strain-induced magnetic field is not perfectly uniform in 
space), whereas in the other cone, the usual continuum 
spectrum would be observed. This motivated us to an- 
alyze the trajectories of a wave packet in such a system 
within the TB model, where we do not need to include 
the pseudo-magnetic fields artificially in the Dirac cones, 
since they appear naturally when we consider the effect 
of the strain-induced changes of the inter-site distances 
on the hopping energies, as explained in the previous sec- 
tion. 

In this subsection, we investigate the dynamics of a 
wave packet with width d ~ 200 A and initial wave vec- 
tor = and fc" = 0.02 A around the Dirac points K 
and K' of Eq. (|2ip in the presence of external and strain- 
induced magnetic field barrier steps. As in the K' valley 
the pseudo-spinor [1, 1]-^ is polarized in the negative y- 
direction of the graphene lattice, we choose [1,-1]^ for 
this case, so that a wave packet in this valley will also 
propagate in the positive ^/-direction. In order to obtain a 
pseudo-magnetic field barrier step, we consider that the 
graphene layer is strained only in the y > region, as 
sketched in Fig. [Sib). We also consider an external mag- 
netic field B — BQ{y)z, where Q{y) is the Heaviside step 
function, which leads to a magnetic barrier step for y > 0, 
described by the vector potential A = {— By <d{y), 0,0). 
In order to avoid effects due to zitterbewegung in the 
(pseudo) magnetic field region, the wave packet starts at 
the position xq — 0, yo = —420 A , so that it can travel 
for some time in the magnetic field-free region y < until 
its velocity converge to a time independent value. 

The infiuence of the external and strain-induced mag- 
netic field barriers on the trajectories of the wave packet 
are analyzed separately in Fig. [71 which shows the trajec- 
tory of the centroid of the wave packets in K (symbols) 
and K' (curves) points, calculated as (r) — {{x), (j/)), (a) 
in a non-strained graphene sheet with magnetic field bar- 
riers B = 5 T (solid, circles) , 7 T (dashed, triangles) and 
10 T (dotted, squares) and (b) in a strained graphene 
sheet with radius R = I fim (solid, circles), 0.8 fim 
(dashed, triangles) and 0.6 fim (dotted, squares). All the 
trajectories form semi-circles in the y > region, which 
is due to the Lorentz force produced by the (pseudo) 
magnetic field. As the external magnetic field (radius 
of the strained region) increases (decreases), the radii 
of these semi-circular trajectories are reduced, since a 
higher (pseudo) magnetic field produces a stronger mod- 
ulus of the Lorentz force. Notice that the radii of tra- 
jectories in the external and pseudo-magnetic fields cases 
are comparable, which means that for radii R = 1 fim - 
0.6 fim of the strained graphene, the generated pseudo- 
magnetic field is also within w 5 T and 10 T. Indeed, the 
strain induced pseudo-magnetic field distribution for the 
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C. Strain induced valley filter 
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FIG. 7: (Color online) Trajectories of the wave packet in the 
x — y plane, obtained by the TB method for such a system, for 
initial momentum ky = 0.02 A ~^ around K (symbols) and K' 
(curves) points, for (a) non-strained graphene with magnetic 
barrier height B = 5 T (solid, circles), 7 T (dashed, triangles) 
and 10 T (dotted, squares), and for (b) a graphene sheet bent 
into an arc of circle with radius R — 1 /im (solid, circles), 0.8 
fim (dashed, triangles) and 0.6 fim (dotted, squares), consid- 
ering _B = T. In (b), symbols and curves coincide for each 
value of R. 
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where /3 w 2 and c is a dimensionless constant which 
depends on the details of the atomic displacements.^^ 
Considering L/i? — >■ in Eq. ([^ the pseudo- magnetic 
field can be approximated as Bs ~ —cf3^o/aR = u/R. 
Using the value w « 4.5 x 10* tA estimated numerically 
in Ref. [l^, one obtains pseudo-magnetic fields within 
Bs « 4.5 T - 7.5 T for ii = 1 ^m - 0.6 /^m, which are 
of the same order of magnitude as the external magnetic 
fields that we considered. For the external magnetic field 
barrier, the trajectories of wave packets in K and K' 
points form circles in opposite directions, as shown in 
Fig. [Tja), which is reasonable, since these packets have 
opposite momentum, which causes a sign change in the 
Lorentz force. Conversely, considering the strain-induced 
magnetic barrier illustrated in Fig. [Hb), the trajectories 
of wave packets in K and K' curve in the same direction, 
since, although their momenta have opposite signs, the 
pseudo-magnetic fields also point in opposite directions 
at each Dirac cone K and K' . 



Let us consider the strained sample in Fig. [6l^b) with 
_R = 1 /im. By comparing the radius of the semi-circular 
trajectory of the wave packet in such a system with those 
obtained for different intensities of the external magnetic 
field barrier, one obtains the strain-induced magnetic 
field for this value of i? as « 4.9 T. Figure El^a) shows 
the trajectories in the x — y plane of the centroid of the 
wave packets in a system where we combine a iZ = 1 /im 
strain for y > Q with an external magnetic field barrier 
S = T (solid, open) and 4.9 T (dashed, full), for wave 
packets in the K (symbols) and K' (curves) Dirac points. 
In the absence of the external magnetic field, both the 
K and K' packets exhibit the same semi-circular trajec- 
tory, as discussed earlier. However, when we combine 
the effect of the strain-induced and external magnetic 
field barriers, the wave packet in K' undergoes a stronger 
Lorentz force and is readily reflected, whereas the one in 
the K point performs a practically straight trajectory, 
as if this packet is not influenced by any Lorentz force. 
This is a consequence of the fact that combining the ef- 
fects of a pseudo-magnetic fleld produced by a i? = l/im 
strain and a B — 4.9 T external magnetic field produces 
a stronger magnetic field in the K' point, while in the 
K point these fields equilibrate, producing a practically 
magnetic field-free region for particles in this cone. In 
this situation, the system works as a valley filter, where 
only wave packets in the K Dirac cone are allowed to pass 
through the strained region, whereas the wave packets in 
K' are reflected. The results for the wave packet in K 
for two other values of the external magnetic field are 
shown as thin solid lines, showing that within a range of 
AB = ±0.2 T around B = 4.9 T, which is a reasonable 
range for magnetic field intensities in experiments, only a 
weak Lorentz force is observed and the valley filter works 
fine. 

The results of Fig. [8] are obtained for both external 
and pseudo-magnetic field barriers starting at the same 
position y = 0. It is straightforward to verify that if 
there is a mismatch between the starting points of the 
strained and external field regions, some deviations will 
occur in the trajectories of the wave packets but, pro- 
vided the length of the mismatch is much smaller than 
the magnetic length, the filtering effect is still stable. As 
an example, a 30 A mismatch between the external and 
pseudo-magnetic field barriers in the system analyzed in 
Fig. [8] would produce a « 5° deviation in the otherwise 
vertical trajectory of the wave packet in K, whereas the 
wave packet in K' is still readily reflected by the combi- 
nation of magnetic fields in the filter region. 

The probability of finding the particle in the 
strained y > region, calculated as 



(23) 



where represents the lines of atomic sites with y > 0, 
is shown as a function of time in Fig. [5^b). In the B = 
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FIG. 8: (Color online) (a) Trajectories on the x — y plane for 
wave packets with initial momentum fc^ = 0.02 around 
K (symbols) and K' (curves) points, considering a graphene 
sheet bent into an arc of circle with radius R = 1 /xm and an 
external magnetic field _B = T (open, solid) and 4.9 T (full, 
dashed). The thin solid curves show the results for two other 
magnetic field intensities for the K packet, (b) Probability of 
finding the particle in y > as a function of time, for wave 
packets with the same configurations as in (a). 



T case, both wave packets in K (open circles) and K' 
(solid) stay in the strained y > region until t « 300 
fs, when they turn back into the y < region, reflected 
by the Lorentz force induced by the strain. However, for 
B = 4.9 T, P> already approaches zero at i « 175 fs for 
the packet in K' (dashed), whereas for K (full circles), it 
remains close to 1 even for large t. 

The efficiency of the proposed valley filter is confirmed 
by Fig. [SI where we present as a function of time for 
initial wave packets given by a combination of Gaussians 
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FIG. 9: (Color online) Probability of finding the electron in 
the filter region y > Q, for an initial wave packet given by a 
combination of Gaussian distributions around both K and K' 
Dirac points, for three different values of the A'-component of 
the wave packet /3. 



around the K and K' points in Eq. ([21 



where k{K') is the Gaussian wave packet in Eq. 
with momentum q around the K{K') Dirac point. The 
results are presented for three difi'erent values of (5, where 
one can easily see that the probability of finding the 
packet in the strained region exhibits a peak at t « 80 fs 
but, as the K' part of the packet is reflected by the mag- 
netic barrier, this probability decays, reaching exactly 
P> — (3 for large t. Such a system proves to be a perfect 
valley fllter, which is able to reflect all the components of 
the incoming packet that are in the K' point and trans- 
mit a wave packet that is composed only of particles with 
momentum in the vicinity of K . 

We point out that when a wave packet reaches the 
edges of a graphene nanoribbon, it can be scattered to a 
different Dirac cone4i Consequently, the efficiency of the 
valley fllter could be compromised if one considers a thin 
nanoribbon, so that the flltered wave packet could still 
reach its boundaries and scatter back to the other valley. 
In order to avoid such an effect, we have considered a 
wide nanoribbon, so that for the time intervals we study 
in this work, boundary effects are not signiflcant. 



D. External and pseudo magnetic field effects on 
the zitterbewegung 

In a previous paper, Rusin and Zawadzki^ used the 
Dirac Hamiltonian for graphene to show that the zitter- 
bewegung, which is transient for P = 0, as discussed ear- 
lier, is permanent for P ^ 0. Furthermore, the authors 
showed that for a Gaussian wave packet, the time evolu- 
tion of the average value of the current is different in x 
and y directions, which they explain as due to the fact 
that the Dirac Hamiltonian is not symmetric in the mo- 
menta Px and Py . Although the same authors say in their 
subsequent paper— that this effect is unphysical, because 
it violates the rotational symmetry of the x — y graphene 
plane, we believe this result is still physical: one must re- 
member that the Dirac model of graphene comes from the 
tight-binding approach for this system, which describes 
a honeycomb lattice of atoms that is not symmetric in 
the x — y plane by deflnition, exhibiting Cg^ symmetry, 
as mentioned in previous section. For each K and K' 
point, the coordinates x and y in the Dirac Hamiltonian 
represent different directions in the real honeycomb lat- 
tice of graphene, where for an infinite sample the x (y) 
coordinate in the Dirac equation is related to one of the 
zigzag (armchair) directions of the real sample. In this 
subsection, we use our TB model of graphene to extend 
the previous study of Rusin and Zawadzki^ to different 
situations. 

We now study the dynamics of a wave packet with 
width d = 200 A , pseudo-spinor ci — 1 and C2 — 1 and 
initial wave vector fc° = = 0, i.e. exactly at one of 
the Dirac points K and K' in Eq. (PT|) . in the presence 
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FIG. 10: (Color online) Electromagnetic dipole radiation as 
a function of time for wave packets with initial pseudo-spinor 
[1, 1]"^ at K (thin curves) and K' (thick curves), considering 
a graphene sheet (a) in the absence of strain and magnetic 
fields, (b) under an uniformly applied magnetic field B = 4.9 
T, (c) bent into an arc of circle with radius R = l/xm (see 
Fig. [G]) and (d) with both the uniform magnetic field B = 4.9 
T and the R — 1/^m bending. Solid (dashed) curves are the 



results for the Sy (ex) component. 



of an uniform applied external magnetic field B = Bz, 
instead of the magnetic field barrier step considered in 
the previous subsection. We also consider a pseudo- 
magnetic field obtained by bending the whole rectangu- 
lar graphene sample into an arc of a circle with radius 
R, as illustrated in Fig. [D^a). The radius of the circle 
is considered as i? = 1 /im, which produces a ~ 4.9 T 
pseudo-magnetic field, as demonstrated in the previous 
subsection. Accordingly, we consider the external uni- 
form magnetic field as i? = 4.9 T. 

A few experimental techniques have been suggested in 
the literature for the observation of zitterbewegung.-iii 
An interesting one^l is based on the fact that the 
wave packet 5'(r, t) exhibits an electric dipole moment 
D{t) = {'i/{r,t)\r\'^{'r,t)) and, consequently, the zitter- 
bewegung yields an oscillation of this dipole moment, 
which is a source of electromagnetic radiation, described 
classically^ by the equation 



m = 



(PD{t) sin$ 



(25) 



where s and $ are the distance and angle of observation, 
respectively, and Cq is the vacuum permittivity. 

Figure [TU] shows the effects of external and strain- 
induced magnetic fields on the electric field radiation 



produced by the zitterbewegung, written in units of 
£c = sin'I>/47reos. Only weak oscillations are expected in 
the armchair (x) direction, since the pseudo-spin [1,1]"^ 
points in the x-direction in the Dirac model which, as 
mentioned earlier, is related to the zigzag (y) direction 
of the honeycomb lattice. Indeed, in Fig. [TUT a-d). the x- 
component of the electric field (dashed) is always close to 
zero. In Fig. [TUT a). we present the results in the absence 
of strain and magnetic fields, for comparison. In this 
case, the oscillations on the electric field are suppressed 
for larger time, which is due to the transient character of 
the zitterbewegung. Besides, the results for Sy (solid) in 
the K (thin curves) and K' (thick curves) points have op- 
posite signs, since for these points the axis of the Dirac 
cones are rotated by angles which differ by 7r/2 differ- 
ence. In Fig. llOf b). the uniformly applied magnetic field 
£? = 4.9 T in an unstrained sample leads to persistent 
oscillations in Sy, which is related to the discrete Lan- 
dau level spectrum created by this field. Each Landau 
level that is populated by the Gaussian distribution con- 
tributes with a different frequency.^ Figure llOf c) shows 
that such a persistent behavior is also obtained by bend- 
ing the graphene sheet into an arc of circle with radius 
R = l^m, which produces a pseudo-magnetic field of 
the same order of magnitude. Notice that the ampli- 
tude of oscillations in this case is four times smaller than 
those found in Fig. [TUTb) for the unstrained sample under 
an external magnetic field. In fact, these two cases are 
not expected to produce the same zitterbewegung, be- 
cause, although both samples exhibit approximately the 
same magnitude of magnetic field, the strained sample 
has not only the pseudo-magnetic field, but also a dif- 
ferent distribution of atomic sites. Thus, in the strained 
case, there is an additional change in the direction of the 
pseudo-spin polarization as the wave packet drifts, due 
to the lattice distortion itself. As we have demonstrated 
in Sec. Ill A, the zitterbewegung strongly depends on 
the pseudo-spin polarization and hence, the different in- 
terplay between the strained atomic sites and the initial 
pseudo-spin polarization produces a different zitterbewe- 
gung for the strained case, as compared to the one of the 
unstrained sample under an external field. 

In Fig. [TUf d). we combine the effects of the R = 1 ^m 
strain and i? = 4.9 T external field to produce a sys- 
tem where the magnetic field is practically zero in the K 
point, but is non-zero in the K' point, so that only the 
packet in the K' point exhibits persistent oscillations. 
For the K point, the external field compensates only the 
effect of the pseudo-magnetic field, namely, the persis- 
tent zitterbewegung, but it does not remove the effect of 
the lattice distortion. As a result, comparing the results 
for K (thin curves) in Figs. [TUT a) and (d), one observes 
that the oscillations are transient in both cases, since 
there is no effective magnetic field, but they still exhibit 
a different oscillation profile at small t, due to the lat- 
tice distortion in the latter case, which is absent in the 
former. 
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IV. CONCLUSION 

We presented a study of the dynamics of Gaussian 
wave packets in graphene under external and strain- 
induced magnetic fields, where the latter is obtained by 
bending the graphene sheet into an arc of a circle. The 
dependence of the zitterbewegung on the initial pseudo- 
spin of the wave packet is investigated, and the results ob- 
tained by means of the tight-binding model and the Dirac 
equation are compared. We demonstrate that the com- 
bination of the pseudo-magnetic field, induced by bend- 
ing the graphene sheet, along with an external magnetic 
field with appropriate strength can be used as an effi- 
cient valley filter. An incoming wave packet composed of 
momenta around the K and K' Dirac points is scattered 
such that all its components in one of the Dirac cones un- 
dergoes a strong Lorentz force and are readily reflected, 
while the components in the other cone are allowed to 
pass through the device with only small distortions in 
their trajectory, due to the very weak residual Lorentz 
force. 

Our results also show that in the absence of exter- 
nal or strain-induced magnetic fields, the zitterbewegung 
is a transient effect, whereas in the presence of any of 



these fields, the oscillations persist in time. In a strained 
sample under an external magnetic field with the appro- 
priate strength, the effective magnetic field in one of the 
Dirac cones is enhanced, whereas in the other cone it 
is practically cancelled. In this situation, a permanent 
zitterbewegung is observed only for wave packets in one 
of the Dirac cones. The wave packet oscillations produce 
electric field radiation, which can be detected experimen- 
tally. 

Finally, we believe the present work contributes to a 
better understanding of the relation between the results 
obtained from TB and Dirac approaches for graphene and 
those to be observed in future experiments on strained 
graphene-based structures. 
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